In [1]:
using Revise
using ShadyPolytopes
using Polyhedra
using Printf
using LaTeXStrings
using SCS
using Latexify
In [2]:
using WGLMakie, Bonito
Page(offline=true)
WGLMakie.activate!()
Inscribed and circumscribed sphere of $\mathcal{C}$¶
In [3]:
function euclidean_norm(x)
return sqrt(sum(t^2 for t in x))
end
function plot_inner_and_outer_sphere(cscb; color=:cornflowerblue)
vertices = all_vertices(cscb)
r = minimum(1/euclidean_norm(w) for w in cscb.normals)
R = maximum(euclidean_norm(v) for v in vertices)
display(@sprintf "R = %.3f" R)
display(@sprintf "r = %.3f" r)
display(@sprintf "R/r = %.3f" R/r)
display("spectral norm bound squared: $(determine_squared_spectralnorm_bound(cscb))")
display(@sprintf "spectral norm bound = %.3f" sqrt(determine_squared_spectralnorm_bound(cscb)))
p = polyhedron(convexhull(convert(Vector{Vector{Float64}},vertices)...))
n = Polyhedra.Mesh(p)
fig,ax,P = mesh(n, transparency=true, alpha=0.7, color=color);
wireframe!(ax, n, color=:black, transparency=true, alpha=0.3, linewidth=1);
outer_sphere = mesh!(ax,Sphere(Point3f(0),R), transparency=true, alpha=0.2, color=:grey)
inner_sphere = mesh!(ax,Sphere(Point3f(0),r), transparency=true, alpha=0.3, color=:black)
return fig
end
display(L"Analyzing $\mathcal{I}$")
display(ShadyPolytopes.perturbed_icosahedron(-3//5,-1//5,1//10))
fig = plot_inner_and_outer_sphere(optimal_icosahedron);
Analyzing $\mathcal{I}$
CSCP{Rational{BigInt}}(Vector{Rational{BigInt}}[[1, -3//5, 1//10], [1, -1//5, 1//10], [1//10, 1, -3//5], [1//10, 1, -1//5], [-3//5, 1//10, 1], [-1//5, 1//10, 1]], Vector{Rational{BigInt}}[[390//1069, -1250//1069, -710//1069], [710//1069, -390//1069, 1250//1069], [20//53, -55//53, 0], [-710//1069, 390//1069, -1250//1069], [-15//17, 0, -20//17], [-10//9, -10//9, -10//9], [-55//53, 0, 20//53], [-1250//1069, -710//1069, 390//1069], [0, -20//53, 55//53], [-20//17, -15//17, 0], [0, 20//17, 15//17], [-390//1069, 1250//1069, 710//1069], [-20//53, 55//53, 0], [0, -20//17, -15//17], [0, 20//53, -55//53], [1250//1069, 710//1069, -390//1069], [10//9, 10//9, 10//9], [20//17, 15//17, 0], [15//17, 0, 20//17], [55//53, 0, -20//53]])
"R = 1.170"
"r = 0.520"
"R/r = 2.253"
"spectral norm bound squared: 137//27"
"spectral norm bound = 2.253"
In [4]:
fig
Out[4]:
Transformation into approximate John position¶
In [5]:
M = approximate_john_position(optimal_icosahedron)
display(latexify(M))
-------------------------------------------------------------
Clarabel.jl v0.10.0 - Clever Acronym
(c) Paul Goulart
University of Oxford, 2022
-------------------------------------------------------------
chordal decomposition:
compact format = on, dual completion = oncronym
(c) Paul Goulart
University of Oxford, 2022
-------------------------------------------------------------
chordal decomposition:
merge method = clique_graph
PSD cones initial = 3
PSD cones decomposable = 1
PSD cones after decomposition = 5
PSD cones after merges = 5
problem:
variables = 22
constraints = 62
nnz(P) = 0
nnz(A) = 238
cones (total) = 9
: Nonnegative = 2, numel = (1,20)
: Power = 2, numel = (3,3)
: PSDTriangle = 5, numel = (6,3,6,10,10)
settings:
linear algebra: direct / qdldl, precision: Float64
max iter = 200, time limit = Inf, max step = 0.990
tol_feas = 1.0e-08, tol_gap_abs = 1.0e-08, tol_gap_rel = 1.0e-08,
static reg : on, ϵ1 = 1.0e-08, ϵ2 = 4.9e-32ym
(c) Paul Goulart
University of Oxford, 2022
-------------------------------------------------------------
chordal decomposition:
merge method = clique_graph
PSD cones initial = 3
PSD cones decomposable = 1
PSD cones after decomposition = 5
PSD cones after merges = 5
problem:
variables = 22
constraints = 62
nnz(P) = 0
nnz(A) = 238
cones (total) = 9
: Nonnegative = 2, numel = (1,20)
: Power = 2, numel = (3,3)
: PSDTriangle = 5, numel = (6,3,6,10,10)
settings:
linear algebra: direct / qdldl, precision: Float64
max iter = 200, time limit = Inf, max step = 0.990
tol_feas = 1.0e-08, tol_gap_abs = 1.0e-08, tol_gap_rel = 1.0e-08,
dynamic reg: on, ϵ = 1.0e-13, δ = 2.0e-07
iter refine: on, reltol = 1.0e-13, abstol = 1.0e-12,
(c) Paul Goulart
University of Oxford, 2022
-------------------------------------------------------------
chordal decomposition:
merge method = clique_graph
PSD cones initial = 3
PSD cones decomposable = 1
PSD cones after decomposition = 5
PSD cones after merges = 5
problem:
variables = 22
constraints = 62
nnz(P) = 0
nnz(A) = 238
cones (total) = 9
: Nonnegative = 2, numel = (1,20)
: Power = 2, numel = (3,3)
: PSDTriangle = 5, numel = (6,3,6,10,10)
settings:
linear algebra: direct / qdldl, precision: Float64
max iter = 200, time limit = Inf, max step = 0.990
tol_feas = 1.0e-08, tol_gap_abs = 1.0e-08, tol_gap_rel = 1.0e-08,
dynamic reg: on, ϵ = 1.0e-13, δ = 2.0e-07
max iter = 10, stop ratio = 5.0
equilibrate: on, min_scale = 1.0e-04, max_scale = 1.0e+04
max iter = 10
iter pcost dcost gap pres dres k/t μ step
---------------------------------------------------------------------------------------------
0 0.0000e+00 -2.4334e+00 2.43e+00 9.34e-01 4.65e+00 1.00e+00 1.00e+00 ------ l Goulart
University of Oxford, 2022
-------------------------------------------------------------
chordal decomposition:
merge method = clique_graph
PSD cones initial = 3
PSD cones decomposable = 1
PSD cones after decomposition = 5
PSD cones after merges = 5
problem:
variables = 22
constraints = 62
nnz(P) = 0
nnz(A) = 238
cones (total) = 9
: Nonnegative = 2, numel = (1,20)
: Power = 2, numel = (3,3)
: PSDTriangle = 5, numel = (6,3,6,10,10)
settings:
linear algebra: direct / qdldl, precision: Float64
max iter = 200, time limit = Inf, max step = 0.990
tol_feas = 1.0e-08, tol_gap_abs = 1.0e-08, tol_gap_rel = 1.0e-08,
dynamic reg: on, ϵ = 1.0e-13, δ = 2.0e-07
max iter = 10, stop ratio = 5.0
equilibrate: on, min_scale = 1.0e-04, max_scale = 1.0e+04
max iter = 10
iter pcost dcost gap pres dres k/t μ step
---------------------------------------------------------------------------------------------
0
1 -2.2412e-01 -5.5373e-01 3.30e-01 1.30e-01 9.48e-01 5.21e-02 1.59e-01 8.61e-01 l Goulart
University of Oxford, 2022
-------------------------------------------------------------
chordal decomposition:
merge method = clique_graph
PSD cones initial = 3
PSD cones decomposable = 1
PSD cones after decomposition = 5
PSD cones after merges = 5
problem:
variables = 22
constraints = 62
nnz(P) = 0
nnz(A) = 238
cones (total) = 9
: Nonnegative = 2, numel = (1,20)
: Power = 2, numel = (3,3)
: PSDTriangle = 5, numel = (6,3,6,10,10)
settings:
linear algebra: direct / qdldl, precision: Float64
max iter = 200, time limit = Inf, max step = 0.990
tol_feas = 1.0e-08, tol_gap_abs = 1.0e-08, tol_gap_rel = 1.0e-08,
dynamic reg: on, ϵ = 1.0e-13, δ = 2.0e-07
max iter = 10, stop ratio = 5.0
equilibrate: on, min_scale = 1.0e-04, max_scale = 1.0e+04
max iter = 10
iter pcost dcost gap pres dres k/t μ step
---------------------------------------------------------------------------------------------
0
1
2 -2.6201e-01 -3.7724e-01 1.15e-01 5.11e-02 3.85e-01 2.00e-02 6.35e-02 6.34e-01
3 -4.1895e-01 -4.4690e-01 2.79e-02 1.49e-02 8.79e-02 4.82e-03 1.61e-02 7.92e-01
4 -4.8033e-01 -4.8614e-01 5.81e-03 3.53e-03 1.76e-02 9.74e-04 3.54e-03 7.92e-01
5 -4.8919e-01 -4.8948e-01 2.86e-04 1.76e-04 8.66e-04 5.00e-05 1.77e-04 9.59e-01
6 -4.8953e-01 -4.8954e-01 7.90e-06 4.89e-06 2.40e-05 1.42e-06 4.92e-06 9.80e-01
7 -4.8954e-01 -4.8954e-01 2.18e-07 1.36e-07 6.66e-07 4.00e-08 1.36e-07 9.80e-01
8 -4.8954e-01 -4.8954e-01 6.05e-09 3.76e-09 1.85e-08 1.12e-09 3.78e-09 9.80e-01
9 -4.8954e-01 -4.8954e-01 1.30e-09 8.08e-10 3.96e-09 2.41e-10 8.12e-10 7.92e-01
---------------------------------------------------------------------------------------------
Terminated with status = solved
solve time = 8.11s
\begin{equation}
\left[
\begin{array}{ccc}
\frac{11}{10} & \frac{11}{10} & \frac{11}{10} \\
0 & \frac{-4}{5} & \frac{4}{5} \\
1 & \frac{-1}{2} & \frac{-2}{5} \\
\end{array}
\right]
\end{equation}
In [6]:
display(L"Analyzing $\mathcal{J} = M \mathcal{I}$")
display([M*x for x in optimal_icosahedron.positive_vertices])
display(optimal_icosahedron_john_position.positive_vertices)
fig = plot_inner_and_outer_sphere(optimal_icosahedron_john_position; color=:yellowgreen)
display(fig)
display("1.55 > spectral_norm_bound: ")
display((155//100)^2 > determine_squared_spectralnorm_bound(optimal_icosahedron_john_position))
Analyzing $\mathcal{J} = M \mathcal{I}$
6-element Vector{Vector{Rational{BigInt}}}:
[11//20, 14//25, 63//50]
[99//100, 6//25, 53//50]
[11//20, -32//25, -4//25]
[99//100, -24//25, -8//25]
[11//20, 18//25, -21//20]
[99//100, 18//25, -13//20]
6-element Vector{Vector{Rational{BigInt}}}:
[11//20, 49//100, 117//100]
[99//100, 57//100, 81//100]
[11//20, 81//100, -51//50]
[99//100, 9//20, -9//10]
[11//20, -13//10, -3//20]
[99//100, -51//50, 9//100]
"R = 1.424"
"r = 0.923"
"R/r = 1.543"
"spectral norm bound squared: 370280747942//155558341125"
"spectral norm bound = 1.543"
"1.55 > spectral_norm_bound: "
true
Determining a projection with small norm¶
In [17]:
function makie_polyhedron(vertices)
vertices = convert(Vector{Vector{Float64}},vertices)
p = polyhedron(convexhull(vertices...))
return Polyhedra.Mesh(p)
end
function plot_projection(cscb; color=:cornflowerblue)
vertices = all_vertices(cscb)
n = makie_polyhedron(vertices)
fig,ax,P = mesh(n, transparency=true, alpha=0.7, color=color);
wireframe!(ax, n, color=:black, transparency=true, alpha=0.3, linewidth=1);
proj, shadiness_const = ShadyPolytopes.find_shadiness_constant_numerically(cscb, false)
approx_proj = ShadyPolytopes.approximate_projection(proj,prec=10^4)
image = makie_polyhedron([1.2*approx_proj*v for v in vertices])
mesh!(ax, image, transparency=true, alpha=0.9, color=:darkred);
kernel = proj*[1,1,0]-[1,1,0]
lines!(ax, [Point3f(kernel),Point3f(-kernel)], transparency=true, alpha=0.5,color=:darkred);
display(@sprintf "s_2 ≈ %.6f" shadiness_const)
display(@sprintf "s_2 ≤ %.6f" ShadyPolytopes.operator_norm(cscb, approx_proj))
display("Projection: ")
display(latexify(approx_proj))
op_norm = ShadyPolytopes.operator_norm(cscb, approx_proj)
display("Operator norm of projection = $op_norm ≈ "*(@sprintf "%.5f" op_norm))
return fig, approx_proj
end
fig, approx_proj = plot_projection(optimal_icosahedron)
fig
presolving:
(round 1, fast) 0 del vars, 1 del conss, 0 add conss, 1 chg bounds, 0 chg sides, 0 chg coeffs, 1 upgd conss, 0 impls, 0 clqs
(0.0s) symmetry computation started: requiring (bin +, int +, cont +), (fixed: bin -, int -, cont -)
(0.0s) symmetry computation finished: 1 generators found (max: 1500, log10 of symmetry group size: 0.0) (symcode time: 0.00)
dynamic symmetry handling statistics:
orbitopal reduction: no components
orbital reduction: no components
lexicographic reduction: no permutations
handled 1 out of 1 symmetry components
(round 2, exhaustive) 0 del vars, 1 del conss, 2 add conss, 1 chg bounds, 0 chg sides, 0 chg coeffs, 1 upgd conss, 0 impls, 0 clqs
presolving (3 rounds: 3 fast, 2 medium, 2 exhaustive):
0 deleted vars, 1 deleted constraints, 0 added constraints, 1 tightened bounds, 0 added holes, 0 changed sides, 0 changed coefficients
0 implications, 0 cliques
presolved problem has 10 variables (0 bin, 0 int, 0 impl, 10 cont) and 132 constraints
123 constraints of type <linear>
9 constraints of type <nonlinear>
Presolving Time: 0.00
time | node | left |LP iter|LP it/n|mem/heur|mdpt |vars |cons |rows |cuts |sepa|confs|strbr| dualbound | primalbound | gap | compl.
0.0s| 1 | 0 | 8 | - | 1826k | 0 | 43 | 132 | 136 | 0 | 0 | 0 | 0 | 1.000000e+00 | -- | Inf | unknown
L 0.2s| 1 | 0 | 8 | - | subnlp| 0 | 43 | 132 | 136 | 0 | 0 | 0 | 0 | 1.000000e+00 | 1.281918e+00 | 28.19%| unknown
0.2s| 1 | 0 | 9 | - | 1834k | 0 | 43 | 132 | 137 | 1 | 1 | 0 | 0 | 1.000000e+00 | 1.281918e+00 | 28.19%| unknown
0.2s| 1 | 0 | 9 | - | 1834k | 0 | 43 | 132 | 137 | 1 | 1 | 0 | 0 | 1.000000e+00 | 1.281918e+00 | 28.19%| unknown
0.2s| 1 | 0 | 1791 | - | 1842k | 0 | 43 | 131 | 137 | 1 | 2 | 0 | 0 | 1.000000e+00 | 1.281918e+00 | 28.19%| unknown
0.2s| 1 | 0 | 1795 | - | 1858k | 0 | 43 | 131 | 141 | 5 | 3 | 0 | 0 | 1.000000e+00 | 1.281918e+00 | 28.19%| unknown
0.2s| 1 | 0 | 1797 | - | 1858k | 0 | 43 | 131 | 144 | 8 | 4 | 0 | 0 | 1.000000e+00 | 1.281918e+00 | 28.19%| unknown
0.2s| 1 | 0 | 1806 | - | 1858k | 0 | 43 | 131 | 149 | 13 | 5 | 0 | 0 | 1.000000e+00 | 1.281918e+00 | 28.19%| unknown
0.2s| 1 | 0 | 1809 | - | 1858k | 0 | 43 | 131 | 152 | 16 | 6 | 0 | 0 | 1.000000e+00 | 1.281918e+00 | 28.19%| unknown
0.2s| 1 | 0 | 1811 | - | 1867k | 0 | 43 | 131 | 154 | 18 | 7 | 0 | 0 | 1.000000e+00 | 1.281918e+00 | 28.19%| unknown
0.2s| 1 | 0 | 1814 | - | 1867k | 0 | 43 | 131 | 156 | 20 | 8 | 0 | 0 | 1.000000e+00 | 1.281918e+00 | 28.19%| unknown
L 0.6s| 1 | 0 | 1814 | - | subnlp| 0 | 43 | 131 | 156 | 20 | 9 | 0 | 0 | 1.000000e+00 | 1.024936e+00 | 2.49%| unknown
0.6s| 1 | 0 | 1814 | - | 1867k | 0 | 43 | 125 | 156 | 20 | 9 | 0 | 0 | 1.000000e+00 | 1.024936e+00 | 2.49%| unknown
0.6s| 1 | 0 | 1819 | - | 1867k | 0 | 43 | 125 | 159 | 23 | 10 | 0 | 0 | 1.000000e+00 | 1.024936e+00 | 2.49%| unknown
0.6s| 1 | 0 | 1822 | - | 1867k | 0 | 43 | 125 | 157 | 27 | 11 | 0 | 0 | 1.000000e+00 | 1.024936e+00 | 2.49%| unknown
time | node | left |LP iter|LP it/n|mem/heur|mdpt |vars |cons |rows |cuts |sepa|confs|strbr| dualbound | primalbound | gap | compl.
0.6s| 1 | 2 | 1822 | - | 1867k | 0 | 43 | 125 | 157 | 27 | 12 | 0 | 0 | 1.000000e+00 | 1.024936e+00 | 2.49%| unknown
L 0.7s| 17 | 14 | 2002 | 121.4 | subnlp| 7 | 43 | 125 | 181 | 133 | 3 | 0 | 0 | 1.000000e+00 | 1.024575e+00 | 2.46%| 0.57%
1.2s| 100 | 47 | 2956 | 29.3 | 2457k | 15 | 43 | 125 | 160 | 750 | 2 | 0 | 0 | 1.000000e+00 | 1.024575e+00 | 2.46%| 4.72%
1.4s| 200 | 69 | 4304 | 21.3 | 3101k | 15 | 43 | 125 | 231 |1542 | 3 | 0 | 0 | 1.000000e+00 | 1.024575e+00 | 2.46%| 13.34%
1.6s| 300 | 111 | 6334 | 21.0 | 4215k | 15 | 43 | 125 | 216 |2856 | 4 | 0 | 0 | 1.000000e+00 | 1.024575e+00 | 2.46%| 13.86%
1.8s| 400 | 153 | 8135 | 20.2 | 4626k | 17 | 43 | 125 | 254 |4017 | 3 | 0 | 0 | 1.000000e+00 | 1.024575e+00 | 2.46%| 17.58%
2.0s| 500 | 161 | 9308 | 18.5 | 4802k | 17 | 43 | 125 | 207 |4709 | 4 | 0 | 0 | 1.000000e+00 | 1.024575e+00 | 2.46%| 24.34%
2.2s| 600 | 213 | 11123 | 18.5 | 5290k | 17 | 43 | 125 | 186 |5884 | 3 | 0 | 0 | 1.000000e+00 | 1.024575e+00 | 2.46%| 30.98%
2.4s| 700 | 233 | 12179 | 17.3 | 6666k | 17 | 43 | 125 | 204 |6502 | 2 | 0 | 0 | 1.000000e+00 | 1.024575e+00 | 2.46%| 31.99%
2.6s| 800 | 271 | 13811 | 17.2 | 7108k | 21 | 43 | 125 | 199 |7471 | 2 | 0 | 0 | 1.000000e+00 | 1.024575e+00 | 2.46%| 32.82%
2.8s| 900 | 287 | 15173 | 16.8 | 7357k | 21 | 43 | 125 | 188 |8349 | 2 | 0 | 0 | 1.000000e+00 | 1.024575e+00 | 2.46%| 38.05%
d 2.9s| 950 | 159 | 16114 | 16.9 |adaptive| 21 | 43 | 125 | 193 |9030 | 3 | 0 | 0 | 1.000000e+00 | 1.012713e+00 | 1.27%| 39.66%
3.1s| 1000 | 155 | 16651 | 16.6 | 7616k | 21 | 43 | 125 | 193 |9364 | 2 | 0 | 0 | 1.000000e+00 | 1.012713e+00 | 1.27%| 46.47%
3.4s| 1100 | 155 | 17850 | 16.2 | 7623k | 21 | 43 | 125 | 217 | 10k| 1 | 0 | 0 | 1.000000e+00 | 1.012713e+00 | 1.27%| 52.29%
3.6s| 1200 | 149 | 19047 | 15.8 | 7639k | 21 | 43 | 125 | 249 | 10k| 2 | 0 | 0 | 1.000000e+00 | 1.012713e+00 | 1.27%| 64.78%
time | node | left |LP iter|LP it/n|mem/heur|mdpt |vars |cons |rows |cuts |sepa|confs|strbr| dualbound | primalbound | gap | compl.
3.8s| 1300 | 153 | 20578 | 15.8 | 7652k | 21 | 43 | 125 | 291 | 12k| 3 | 0 | 0 | 1.000000e+00 | 1.012713e+00 | 1.27%| 68.96%
3.9s| 1400 | 105 | 21422 | 15.3 | 7662k | 21 | 43 | 125 | 0 | 12k| 0 | 0 | 0 | 1.008217e+00 | 1.012713e+00 | 0.45%| 88.43%
4.1s| 1500 | 29 | 21974 | 14.6 | 7688k | 21 | 43 | 125 | 266 | 13k| 1 | 0 | 0 | 1.011336e+00 | 1.012713e+00 | 0.14%| 98.38%
* 4.2s| 1537 | 6 | 22137 | 14.4 | LP | 24 | 43 | 122 | 350 | 13k| 2 | 0 | 0 | 1.012657e+00 | 1.012713e+00 | 0.01%| 99.99%
SCIP Status : problem is solved [optimal solution found]
Solving Time (sec) : 4.25
Solving Nodes : 1547
Primal Bound : +1.01271284807885e+00 (5 solutions)
Dual Bound : +1.01271284807885e+00
Gap : 0.00 %
"s_2 ≈ 1.012713"
"s_2 ≤ 1.012747"
"Projection: "
\begin{equation}
\left[
\begin{array}{ccc}
\frac{82602121}{79729122} & \frac{54836807}{79729122} & \frac{-722323}{13288187} \\
\frac{-4217717}{79729122} & \frac{-774259}{79729122} & \frac{1060409}{13288187} \\
\frac{695635}{39864561} & \frac{13277555}{39864561} & \frac{12938397}{13288187} \\
\end{array}
\right]
\end{equation}
"Operator norm of projection = 14386149522//14205071903 ≈ 1.01275"
Out[17]:
In [ ]: